Confidence in Education

Author

Ari Gelbard


This document aims to explores the relationship between confidence in the US education system and actual government spending on education.

Loading Packages and Data

In this section, we pull in data downloaded from the Bureau of Economic Analysis (BEA) website and from the General Social Survey (GSS).

### Load packages
library(tidyverse)
library(readxl)
require(forecast)
library(zoo)
library(kableExtra)
library(astsa)

### Load data
gved_y_raw <- read_excel("data/nipa_3155.xlsx")
gv_y_raw   <- read_excel("data/nipa_116a.xlsx")
gv_q_raw   <- read_excel("data/nipa_116q.xlsx")
pop_raw    <- read_excel("data/nipa_21q.xlsx")
gss_raw    <- read_excel("data/gss_educ.xlsx")

Cleaning data

First, lets clean the data on government spending over time in the US. Unfortunately, the BEA’s data on inflation adjusted spending on education only goes back to 2007. To work around this, we can take the ratio of nominal education spending to nominal total government spending and multiply real total government spending to get an estimate of real total government spending on education. However, total government spending on education is still misleading since it is increasing with things like population growth. To fix this, we can calculate education spending per capita by dividing by the U.S. population. We create two dataframes for the analysis: one with yearly data that we can use to compare to the yearly GSS data, and another that is quarterly for time series analysis.

### BEA NIPA data 
## Population to get per capita spending 
pop <- pop_raw %>%
  select(-1) %>%
  slice(-(1:4)) %>%
  filter(row_number() %in% c(1,2,45)) %>%
  t() %>% as.data.frame() %>%
  slice(-1) %>%
  rename(year = V1, quarter = V2, popq = V3) %>%
  fill(year, .direction = "down") %>% 
  mutate(yearq = paste0(year,quarter),
         year  = as.numeric(year),
         popq = as.numeric(popq)) %>%
  select(year, yearq, popq)

pop_y <- pop %>%
  group_by(year) %>%
  summarise(popy = mean(popq, na.rm = TRUE))

## Getting ratio of education spending to total
govt_educ_y <- gved_y_raw %>%
  select(-1) %>%
  slice(-(1:4)) %>%
  filter(row_number() %in% c(1,3,31)) %>%
  t() %>% as.data.frame() %>%
  slice(-1) %>%
  rename(year = V1, govt = V2, educ = V3) %>%
  mutate(educ = as.numeric(educ),
         govt = as.numeric(govt),
         year = as.numeric(year)) %>%
  mutate(educ_govt = educ / govt)

## Getting annual government spending 
govt_y <- gv_y_raw %>%
  select(-1) %>%
  slice(-(1:4)) %>%
  filter(row_number() %in% c(1,24)) %>% 
  t() %>% as.data.frame() %>%
  slice(-1) %>%
  rename(year = V1, govty = V2) %>%
  mutate(year = as.numeric(year),
         govty = as.numeric(govty)) %>%
  left_join(govt_educ_y %>% select(year, educ_govt), by = "year") %>%
  filter(!is.na(educ_govt)) %>%
  mutate(gved_y = as.numeric(govty) * as.numeric(educ_govt)) %>%
  inner_join(pop_y, by = "year") %>%
  mutate(gved_pc_y = (gved_y * 1e6) / (popy * 1e3)) %>%
  select(year, gved_pc_y)
  
## Quarterly government spending on education (for time series forecast)
govt_q <- gv_q_raw %>%
  select(-1) %>%
  slice(-(1:4)) %>%
  filter(row_number() %in% c(1,2,24)) %>% 
  t() %>% as.data.frame() %>%
  slice(-1) %>%
  rename(year = V1, quarter = V2, govtq = V3) %>%
  fill(year, .direction = "down") %>% 
  mutate(yearq   = paste0(year,quarter),
         year    = as.numeric(year),
         quarter = as.numeric(substr(quarter,2,2)),
         govtq   = as.numeric(govtq)) %>%
  select(c(year,quarter,yearq,govtq)) %>%
  left_join(govt_educ_y %>% select(year, educ_govt), by = "year") %>%
  filter(!is.na(educ_govt)) %>%
  mutate(gved_q = govtq * educ_govt) %>%
  inner_join(pop, by = c("yearq","year")) %>%
  mutate(gved_pc_q = (gved_q * 1e6) / (popq * 1e3),
         yearq = as.yearqtr(yearq)) %>%
  select(c(year,quarter,yearq, gved_pc_q))

Regarding the variables from the GSS, we use the following surveys:

  • The coneduc survey asks whether the respondent has a great deal of confidence, only some confidence, or hardly any confidence in the institution of education in the US.
  • The nateduc survey asks if the respondent thinks the US government is spending too much, too little, or the right amount on education.

We create two dataframes, one that has the average response per year and another that has the number of responses per response per year.

### GSS Education variables
## Taking an average of responses per year
gss_num <- gss_raw %>%
  mutate(nateduc_num = case_when(
            nateduc == "TOO LITTLE" ~ 1,
            nateduc == "ABOUT RIGHT" ~ 2,
            nateduc == "TOO MUCH" ~ 3,
            TRUE ~ NA_real_
          ),
         coneduc_num = case_when(
            coneduc == "HARDLY ANY" ~ 1,
            coneduc == "ONLY SOME" ~ 2, 
            coneduc == "A GREAT DEAL" ~ 3,
            TRUE ~ NA_real_
            ),
         year = as.numeric(year)) %>%
  group_by(year) %>% 
  summarise(mean_nateduc = mean(nateduc_num, na.rm = TRUE),
            mean_coneduc = mean(coneduc_num, na.rm = TRUE),
            n_obs = n()) %>%
  ungroup() %>%
  mutate(mean_coneduc = na.approx(mean_coneduc, year, na.rm = FALSE))

## Counting how many responses of each type per year
gss_count <- gss_raw %>%
  mutate(nateduc = ifelse(nateduc %in% c("TOO LITTLE","ABOUT RIGHT","TOO MUCH"), nateduc, NA),
         coneduc = ifelse(coneduc %in% c("A GREAT DEAL","HARDLY ANY","ONLY SOME"), coneduc, NA)) %>%
  select(-id_) %>%
  mutate(year = as.numeric(year))  %>%
  pivot_longer(
    cols = c(nateduc, coneduc),
    names_to = "question",
    values_to = "response"
  ) %>%
  group_by(year, question, response) %>%
  summarise(n = n(), .groups = "drop")

Analysis

First, lets look at how government spending on education per capita has evolved over time.

## NIPA Data
ggplot(data = govt_y, aes(x=year,y=gved_pc_y)) +
  geom_line() + 
  labs(title="Real Government Spending on Education Per Capita",
       x = "Year",
       y = "Chained (2017) Dollars") +
  theme_minimal()

This chart shows that the government has been spending more per capita on education over time. There were some dips in the 1980s and in the early 2010s, but overall it has been increasing. However, spending has stagnated a bit since the early 2000s.

Now that we’ve seen how spending has changed over time, lets look at confidence in education. We can make line charts for the proportion of each response the coneduc question.

## GSS Data
gss_prop <- gss_count %>%
  filter(!is.na(response)) %>%
  group_by(year, question) %>%
  mutate(prop = n / sum(n)) 

# Confidence in Education
coneduc_prop <- gss_prop %>%
  filter(question == "coneduc")

ggplot(coneduc_prop, aes(x = year, y = prop, color = response)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    title = "How much confidence do you have in the education\nsystem in the US?",
    x = "Year",
    y = "Proportion of Respondents",
    color = "Response"
  ) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 0, vjust = 0.5))

It appears that confidence in the education system has decreased over time. The proportion of respondents who have a great deal of confidence has decreased from around 40% to less than 20% and the proportion who have hardly any confidence has increased from less than 10% to around 25%.

We can also perform a chi-squared test for independence to test if the distribution of GSS responses is different between two time periods. The hypotheses are:

\(H_0:\) The distribution of responses to how much confidence in education is the same prior to 1990 vs after 1990.

\(H_a:\) The distribution of responses is not the same prior to 1990 vs after 1990.

# Choose two periods for comparison (example: early vs recent)
gss_test1 <- gss_raw %>%
  filter(coneduc %in% c("HARDLY ANY","ONLY SOME","A GREAT DEAL")) %>%
  mutate(period = ifelse(year <= 1990, "Early", "Late")) %>%
  group_by(period, coneduc) %>%
  summarise(n = n(), .groups = "drop") %>%
  pivot_wider(names_from = coneduc, values_from = n)

gss_test1
# A tibble: 2 × 4
  period `A GREAT DEAL` `HARDLY ANY` `ONLY SOME`
  <chr>           <int>        <int>       <int>
1 Early            7076         2366       11313
2 Late             7206         5488       17115
# Run the chi-square test
chisq.test(gss_test1[ , -1])   

    Pearson's Chi-squared test

data:  gss_test1[, -1]
X-squared = 831.81, df = 2, p-value < 2.2e-16

Since the p-value is less than 0.05, we can reject the null hypothesis at the 5% level. In other words, we have sufficient statistical evidence to conclude that the distribution of responses is different prior to 1990 vs after 1990.

What about feelings about education spending? Let’s make a similar line chart over time to explore this.

# Spending on Education
nateduc_prop <- gss_prop %>%
  filter(question == "nateduc")

ggplot(nateduc_prop, aes(x = year, y = prop, color = response)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    title = "Is US spending on education too much, too little, or about right?",
    x = "Year",
    y = "Proportion of Respondents",
    color = "Response"
  ) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 0, vjust = 0.5))

Since the beginning of the survey, most respondents believed that the government was spending too little on education, and this proportion has increased significantly to almost 80% of respondents.

We can also perform a chi-squared test for independence to test if the distribution of GSS responses is different between two time periods. The hypothese are:

\(H_0:\) The distribution of responses to how much are we spending on education is the same prior to 1990 vs after 1990.

\(H_a:\) The distribution of responses is not the same prior to 1990 vs after 1990.

# Choose two periods for comparison (example: early vs recent)
gss_test2 <- gss_raw %>%
  filter(nateduc %in% c("TOO LITTLE","ABOUT RIGHT","TOO MUCH")) %>%
  mutate(period = ifelse(year <= 1990, "Early", "Late")) %>%
  group_by(period, nateduc) %>%
  summarise(n = n(), .groups = "drop") %>%
  pivot_wider(names_from = nateduc, values_from = n)

gss_test2
# A tibble: 2 × 4
  period `ABOUT RIGHT` `TOO LITTLE` `TOO MUCH`
  <chr>          <int>        <int>      <int>
1 Early           6183        10351       1440
2 Late            4989        17102       1432
# Run the chi-square test
chisq.test(gss_test2[ , -1])   

    Pearson's Chi-squared test

data:  gss_test2[, -1]
X-squared = 1064.8, df = 2, p-value < 2.2e-16

Again, since the p-value is less than 0.05, we can reject the null hypothesis at the 5% level. In other words, we have sufficient statistical evidence to conclude that the distribution of responses is different prior to 1990 vs after 1990.

Finally, lets plot the mean of these two variables over time to create an even more clear visualization of the decline in sentiment.

gss_long <- gss_num %>%
  pivot_longer(
    cols = c(mean_nateduc, mean_coneduc),
    names_to = "variable",
    values_to = "mean_value"
  )

ggplot(gss_long, aes(x = year, y = mean_value, color = variable)) +
  geom_line(size = 1.1) +
  labs(
    title = "Mean Education Attitudes Over Time",
    x = "Year",
    y = "Mean Value",
    color = "Variable"
  ) +
  theme_minimal()

Now that we’ve looked at the variables of interest separately, lets investigate how spending and confidence in education relate to each other. The following charts show how opinions about the amount of education spending relate to actual education spending.

nateduc_nipa <- gss_prop %>%
  filter(question == "nateduc") %>%
  left_join(govt_y, by = "year")

ggplot(nateduc_nipa, aes(x = gved_pc_y, y = prop)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_wrap(~ response) +
  labs(
    title = "Confidence in amount of Education Spending vs Actual Spending",
    x = "Real Education Spending Per Capita",
    y = "Proportion of Respondents"
  ) +
  scale_y_continuous(labels = scales::percent) +
  theme_bw()
`geom_smooth()` using formula = 'y ~ x'

The negative slope for the respondents who believe the government spends about the right amount on education and the positive slope for the those who believe the government spends too little on education indicate that even when the government prioritize education more, people still believe that the government isn’t spending enough. This may indicate that despite some efforts to spend more, the government still isn’t spending as much on education as they should.

Doing the same with confidence in the education system gives similar results.

coneduc_nipa <- gss_prop %>%
  filter(question == "coneduc") %>%
  left_join(govt_y, by = "year")

ggplot(coneduc_nipa, aes(x = gved_pc_y, y = prop)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_wrap(~ response) +
  labs(
    title = "Confidence in the US Education System vs Spending",
    x = "Real Education Spending Per Capita",
    y = "Proportion of Respondents"
  ) +
  scale_y_continuous(labels = scales::percent) +
  theme_bw()
`geom_smooth()` using formula = 'y ~ x'

Finally, we can do the same thing but using the mean of the surveys.

gss_num_govt <- gss_num %>%
  inner_join(govt_y, by = "year") %>%
  slice(-1)

# Mean_nateduc vs govt spending
ggplot(gss_num_govt, aes(x = gved_pc_y, y = mean_nateduc)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(
    title = "National Education Spending Attitudes vs Education Spending",
    x = "Govt Spending per Capita",
    y = "Mean National Education Attitude"
  ) +
  theme_minimal()
`geom_smooth()` using formula = 'y ~ x'

# Mean_coneduc vs govt spending
ggplot(gss_num_govt, aes(x = gved_pc_y, y = mean_coneduc)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(
    title = "Confidence in Education vs Education Spending",
    x = "Govt Spending per Capita",
    y = "Mean Confidence in Education"
  ) +
  theme_minimal()
`geom_smooth()` using formula = 'y ~ x'

All these plots serve to show that even as government spending on education increases, people still increasingly have lost trust in the education system in the US. Let’s explore how these trends may evolve in the future with time series forecasts.

ARIMA

In this section, we analyze and attempt to forecast the time series for education spending and confidence in education. We focus only on confidence in education rather than feelings about speding on education since it more closely relates to our overall topic of trust in american institutions.

Government Spending on Education

gv_ts <- ts(govt_q$gved_pc_q,
            start = c(min(govt_q$year),min(govt_q$quarter)),
            frequency = 4) 

Lets start with a classical decomposition.

gv_ts_decomp <- decompose(gv_ts)
autoplot(gv_ts_decomp)

The trend component of the chart shows a clear upward trend over time that matches the data. The seasonal component is also showing clear seasonality with spikes that appear to occur at the beginning of each year, however the effect of the seasonality seems small compared to the trend. Finally, the remainder shows fluctuations around zero that aren’t very uniform. It doesn’t seem likely that the remainder is showing cyclical patterns like recessions since it doesn’t seem to follow any particular pattern, but perhaps it reflects other outside political pressures on spending that are difficult to predict.

Let’s plot the lag and autocorrelation function plots to investigate the patterns more closely.

gglagplot(gv_ts, do.lines=FALSE, lags = 9) +
  xlab("Lags") + ylab("Yt") + 
  ggtitle("Lag Plot for education spending") + 
  theme(axis.text.x=element_text(angle=45, hjust=1))

The linear relationships in the lag plots suggest high serial correlation among consecutive observations.

ggAcf(gv_ts, 50) +
  ggtitle("Autocorrelation Function for Real\nGovernment Spending on Education Per Capita") +
  theme_minimal()

This plot shows clear autocorrleation in the data, however the seasonality isn’t as present. Lets confirm stationarity with an augmented Dickey Fuller test.

tseries::adf.test(gv_ts)

    Augmented Dickey-Fuller Test

data:  gv_ts
Dickey-Fuller = -2.4123, Lag order = 6, p-value = 0.4023
alternative hypothesis: stationary

Because the p-values for this test is greater than 0.05, we cannot reject the null hypothesis that these series are not stationary, confirming our suspicions from the ACF plot.

We can try to correct this with differencing.

diff_1 <- diff(gv_ts)
ggAcf(diff_1,50) +
  ggtitle("ACF of First-Order Differenced Series") +
  theme_minimal()

It appears that one round of differencing may not have been sufficient since there is still a decreasing pattern every 4 lags. Let’s try one more time to be safe.

diff_2 <- diff(gv_ts, differences = 2)
ggAcf(diff_2,50) +
  ggtitle("ACF of Second-Order Differenced Series") +
  theme_minimal()

ggPacf(diff_2,50) +
  ggtitle("PACF of the Differenced Series") +
  theme_minimal()

This looks a little bit better, perhaps we will need to test for both 1 and 2 lags in the Arima. With this, we can proceed with the ARIMA.

First lets find the optimal parameters. Based on the differenced ACF and PACF plot, I believe the correct range for q is 0-8 and for p is 0-4.

# Define parameter ranges
p_range <- 0:4
d_range <- 1:2
q_range <- 0:9

# Calculate total combinations
n_combinations <- length(p_range) * length(d_range) * length(q_range)

# Create an empty matrix to store results
results_matrix <- matrix(NA, nrow = n_combinations, ncol = 6)

# Initialize index for matrix row
i <- 1

# Loop through combinations of ARIMA model parameters
for (d in d_range) {
  for (q in q_range) {
    for (p in p_range) {
  
      # Fit ARIMA model with specified (p,d,q)
      model <- Arima(gv_ts, order = c(p, d, q), include.drift = TRUE)
  
      # Store model parameters and AIC/BIC/AICc values in matrix
      results_matrix[i, ] <- c(p, d, q, model$aic, model$bic, model$aicc)
  
      # Increment row index
      i <- i + 1
    }
  }
}

# Convert matrix to data frame
results_df <- as.data.frame(results_matrix)
colnames(results_df) <- c("p", "d", "q", "AIC", "BIC", "AICc")

# Find the row with the lowest AIC
highlight_row <- which.min(results_df$AIC)

# Generate kable table with highlighting for the row with the lowest AIC
knitr::kable(results_df, align = 'c', caption = "Comparison of ARIMA Models") %>%
  kable_styling(full_width = FALSE, position = "center") %>%
  row_spec(highlight_row, bold = TRUE, background = "#FFFF99")  # Highlight row in yellow
Comparison of ARIMA Models
p d q AIC BIC AICc
0 1 0 2491.868 2499.012 2491.914
1 1 0 2489.464 2500.181 2489.557
2 1 0 2491.169 2505.458 2491.324
3 1 0 2493.105 2510.965 2493.338
4 1 0 2413.508 2434.941 2413.836
0 1 1 2489.862 2500.578 2489.955
1 1 1 2477.350 2491.639 2477.506
2 1 1 2475.754 2493.614 2475.987
3 1 1 2469.481 2490.914 2469.809
4 1 1 2413.555 2438.560 2413.995
0 1 2 2491.114 2505.402 2491.269
1 1 2 2461.185 2479.046 2461.419
2 1 2 2471.032 2492.465 2471.361
3 1 2 2450.991 2475.996 2451.430
4 1 2 2415.378 2443.955 2415.945
0 1 3 2479.648 2497.509 2479.881
1 1 3 2453.874 2475.306 2454.202
2 1 3 2453.964 2478.969 2454.403
3 1 3 2448.673 2477.250 2449.240
4 1 3 2416.449 2448.599 2417.161
0 1 4 2426.689 2448.122 2427.017
1 1 4 2428.569 2453.574 2429.009
2 1 4 2427.455 2456.032 2428.022
3 1 4 2429.022 2461.172 2429.734
4 1 4 2418.418 2454.140 2419.291
0 1 5 2428.602 2453.607 2429.041
1 1 5 2424.873 2453.450 2425.440
2 1 5 2426.905 2459.054 2427.616
3 1 5 2427.978 2463.700 2428.851
4 1 5 2419.521 2458.815 2420.573
0 1 6 2429.563 2458.141 2430.130
1 1 6 2426.622 2458.771 2427.333
2 1 6 2423.968 2459.690 2424.841
3 1 6 2425.868 2465.161 2426.919
4 1 6 2421.496 2464.362 2422.744
0 1 7 2430.482 2462.632 2431.194
1 1 7 2428.499 2464.220 2429.372
2 1 7 2425.761 2465.055 2426.813
3 1 7 2427.382 2470.248 2428.630
4 1 7 2423.327 2469.765 2424.789
0 1 8 2426.343 2462.065 2427.216
1 1 8 2423.969 2463.263 2425.021
2 1 8 2425.870 2468.736 2427.118
3 1 8 2420.360 2466.798 2421.822
4 1 8 2423.777 2473.787 2425.471
0 1 9 2425.562 2464.856 2426.614
1 1 9 2425.747 2468.613 2426.995
2 1 9 2428.923 2475.361 2430.385
3 1 9 2422.604 2472.615 2424.298
4 1 9 2420.926 2474.508 2422.869
0 2 0 2694.714 2698.283 2694.730
1 2 0 2589.678 2596.815 2589.724
2 2 0 2556.299 2567.004 2556.392
3 2 0 2421.927 2436.201 2422.083
4 2 0 2418.246 2436.088 2418.481
0 2 1 2483.934 2491.071 2483.980
1 2 1 2471.136 2481.841 2471.229
2 2 1 2466.091 2480.365 2466.247
3 2 1 2416.039 2433.881 2416.274
4 2 1 2409.633 2431.043 2409.962
0 2 2 2461.556 2472.261 2461.649
1 2 2 2462.294 2476.568 2462.450
2 2 2 2458.668 2476.510 2458.902
3 2 2 2416.269 2437.680 2416.599
4 2 2 2409.464 2434.442 2409.905
0 2 3 2461.613 2475.886 2461.769
1 2 3 2463.498 2481.340 2463.732
2 2 3 2448.353 2469.763 2448.683
3 2 3 2418.269 2443.248 2418.710
4 2 3 2411.237 2439.784 2411.807
0 2 4 2460.940 2478.781 2461.174
1 2 4 2444.143 2465.554 2444.473
2 2 4 2446.544 2471.523 2446.985
3 2 4 2418.292 2446.839 2418.861
4 2 4 2412.406 2444.521 2413.121
0 2 5 2421.724 2443.134 2422.053
1 2 5 2422.362 2447.340 2422.803
2 2 5 2424.094 2452.640 2424.663
3 2 5 2420.282 2452.397 2420.996
4 2 5 2414.341 2450.025 2415.218
0 2 6 2422.448 2447.427 2422.889
1 2 6 2418.950 2447.497 2419.519
2 2 6 2420.950 2453.065 2421.664
3 2 6 2420.879 2456.562 2421.755
4 2 6 2415.589 2454.841 2416.645
0 2 7 2424.378 2452.924 2424.947
1 2 7 2420.950 2453.065 2421.664
2 2 7 2423.793 2459.477 2424.670
3 2 7 2420.419 2459.671 2421.475
4 2 7 2417.589 2460.409 2418.842
0 2 8 2426.250 2458.365 2426.964
1 2 8 2422.027 2457.710 2422.903
2 2 8 2422.764 2462.015 2423.820
3 2 8 2414.724 2457.545 2415.978
4 2 8 2419.284 2465.673 2420.752
0 2 9 2422.732 2458.416 2423.609
1 2 9 2420.357 2459.609 2421.413
2 2 9 2423.369 2466.189 2424.622
3 2 9 2415.815 2462.204 2417.283
4 2 9 2419.702 2469.659 2421.402
auto.arima(gv_ts)
Series: gv_ts 
ARIMA(1,2,1)(0,0,1)[4] 

Coefficients:
          ar1      ma1    sma1
      -0.2286  -0.9327  0.4306
s.e.   0.0642   0.0314  0.0518

sigma^2 = 581.6:  log likelihood = -1205.54
AIC=2419.09   AICc=2419.24   BIC=2433.36

The parameter search suggests an ARIMA(4,2,2) while the auto.arima() suggests ARIMA(1,2,1). Lets compare these to determine which one would be better.

model_output1 <- capture.output(sarima(gv_ts,4,2,2))

# Find the line numbers dynamically based on a keyword
start_line <- grep("Coefficients", model_output1)  # Locate where coefficient details start
end_line <- length(model_output1)  # Last line of output

# Print the relevant section automatically
cat(model_output1[start_line:end_line], sep = "\n")
Coefficients: 
    Estimate     SE t.value p.value
ar1   0.0045 0.1008  0.0451  0.9641
ar2   0.0385 0.0536  0.7171  0.4740
ar3   0.0810 0.0530  1.5302  0.1272
ar4   0.5197 0.0540  9.6317  0.0000
ma1  -1.1826 0.1191 -9.9301  0.0000
ma2   0.1826 0.1182  1.5449  0.1236

sigma^2 estimated as 536.3868 on 256 degrees of freedom 
 
AIC = 9.196427  AICc = 9.197684  BIC = 9.291764 
 
model_output2 <- capture.output(sarima(gv_ts,1,2,1))

# Find the line numbers dynamically based on a keyword
start_line <- grep("Coefficients", model_output2)  # Locate where coefficient details start
end_line <- length(model_output2)  # Last line of output

# Print the relevant section automatically
cat(model_output2[start_line:end_line], sep = "\n")
Coefficients: 
    Estimate     SE  t.value p.value
ar1  -0.2596 0.0652  -3.9825   1e-04
ma1  -0.8792 0.0360 -24.4325   0e+00

sigma^2 estimated as 708.7791 on 260 degrees of freedom 
 
AIC = 9.431818  AICc = 9.431995  BIC = 9.472677 
 

Both Residual Plots show nearly consistent fluctuation around zero, suggesting that the residuals are nearly stationary with a constant mean and finite variance over time.

Both Autocorrelation Functions (ACF) of the residuals reveal very few significant autocorrelations, however the plot for the ARIMA(1,2,1) does have a few more autocorrelations than the ARIMA(4,2,2).

Both Q-Q Plots indicate that the residuals follow a near-normal distribution, with minor deviations at the tails, which is typical in time series data.

The Ljung-Box Test p-values for the ARIMA(4,2,2) remain above the 0.05 significance level, implying no significant autocorrelations left in the residuals while in the ARIMA(1,2,1), all p-values remain below the 0.05 significance level implying some significant autocorrelations left in the residuals and concluding that model improvements may necessary.

Coefficient Significance: The model coefficients in the ARIMA(1,2,1) model are all significant while the coefficients for the first three AR terms and the last MA term are insignificant in the ARIMA(4,2,2) model.

fit1 <- Arima(gv_ts, order = c(4,2,2), include.drift = FALSE)
summary(fit1)
Series: gv_ts 
ARIMA(4,2,2) 

Coefficients:
         ar1     ar2    ar3     ar4      ma1     ma2
      0.0045  0.0385  0.081  0.5197  -1.1826  0.1826
s.e.  0.1008  0.0536  0.053  0.0540   0.1191  0.1182

sigma^2 = 549:  log likelihood = -1197.73
AIC=2409.46   AICc=2409.9   BIC=2434.44

Training set error measures:
                    ME     RMSE      MAE         MPE      MAPE      MASE
Training set -1.015172 23.07227 17.61682 -0.03459217 0.8088008 0.3449994
                     ACF1
Training set -0.002555783
fit2 <- Arima(gv_ts, order = c(1,2,1), include.drift = FALSE)
summary(fit2)
Series: gv_ts 
ARIMA(1,2,1) 

Coefficients:
          ar1      ma1
      -0.2596  -0.8792
s.e.   0.0652   0.0360

sigma^2 = 714.2:  log likelihood = -1232.57
AIC=2471.14   AICc=2471.23   BIC=2481.84

Training set error measures:
                    ME     RMSE      MAE        MPE      MAPE      MASE
Training set 0.2117617 26.52199 20.79503 0.01694402 0.9670959 0.4072402
                    ACF1
Training set -0.03802343

Looking at the error measurements, both the MAE and the RMSE are smaller in the ARIMA(4,2,2) model.

naive <- naive(gv_ts,h=5)
mean  <- meanf(gv_ts,h=5)
drift <- rwf(gv_ts,drift=TRUE,h=5)
accuracy(naive)
                   ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
Training set 8.835734 28.79736 21.69429 0.4615319 1.017407 0.4248508 -0.1290414
accuracy(mean)
                        ME     RMSE      MAE      MPE     MAPE     MASE
Training set -9.354218e-14 661.1524 594.0957 -11.0776 31.14832 11.63449
                  ACF1
Training set 0.9869592
accuracy(drift)
                        ME     RMSE      MAE        MPE      MAPE      MASE
Training set -1.085001e-13 27.40835 21.14088 0.04062096 0.9861303 0.4140131
                   ACF1
Training set -0.1290414

Furthermore, both models beat the benchmark models. Overall, it seems like the ARIMA(4,2,2) model is a better fit.

Now we can finally generate our forecast with the ARIMA(4,2,2) model.

# Generate the forecast
forecast_result1 <- forecast(fit1, h = 20)

# Plot the forecast
autoplot(forecast_result1) +
  labs(title = "ARIMA(4,2,2) Forecast",
       x = "Time",
       y = "Real Education Spending per Capita") +
  theme_minimal()

This forecast suggests that real education spending per capita will continue to increase over time.

Educational Confidence

Now lets analyze confidence in the education system, since this is our primary variable of interest.

coneduc_ts_df <- gss_num %>%
  select(year, mean_coneduc) %>%
  arrange(year) %>%
  complete(year = full_seq(year, 1)) %>% 
  mutate(mean_coneduc = na.approx(mean_coneduc, year, na.rm = FALSE)) %>%
  slice(-1)

coneduc_ts <- ts(coneduc_ts_df$mean_coneduc,
                 start = min(coneduc_ts_df$year),
                 end   = max(coneduc_ts_df$year),
                 frequency = 1)

Since this data is annual, we can’t use the decompose() function for a classical decomposition. Instead, we will extract the trend and remainder parts separately.

### Trend 
trend <- stats::lowess(coneduc_ts)$y
remainder <- coneduc_ts - trend

# Plots
plot(coneduc_ts, type = "l",
     main = "Data",
     xlab = "Time", ylab = "Confidence in Education")

plot(trend, type = "l",
     main = "Trend Component",
     xlab = "Time", ylab = "Trend")

plot(remainder, type = "l",
     main = "Remainder Component",
     xlab = "Time", ylab = "Remainder")

The trend component shows a clear downward trend in sentiment and the remainder shows some fluctuations around zero with some small spikes in the early years of the survey.

gglagplot(coneduc_ts, do.lines=FALSE, lags = 9) +
  xlab("Lags") + ylab("Yt") + 
  ggtitle("Lag Plot for confidence in education") + 
  theme(axis.text.x=element_text(angle=45, hjust=1))

The lag plots show a slight relationship between consecutive observations, suggesting there may not be very much autocorrelation.

ggAcf(coneduc_ts, 50) + 
  ggtitle("Autocorrelation Function for\nAverage Confidence in the Education System") + 
  theme_minimal()

The ACF plot, however, shows some autocorrleation. Lets confirm stationarity with an augmented Dickey Fuller test.

tseries::adf.test(coneduc_ts)

    Augmented Dickey-Fuller Test

data:  coneduc_ts
Dickey-Fuller = -2.5281, Lag order = 3, p-value = 0.362
alternative hypothesis: stationary

Because the p-values for this test is greater than 0.05, we cannot reject the null hypothesis that these series are not stationary, confirming our conclusions from the ACF plot.

Lets try to correct this with differencing.

diff_1 <- diff(coneduc_ts)
ggAcf(diff_1,50) +
  ggtitle("ACF of First-Order Differenced Series") +
  theme_minimal()

ggPacf(diff_1,50) +
  ggtitle("PACF of First-Order Differenced Series") +
  theme_minimal()

The ACF and PACF plots for the differenced series indicate that the data is now stationary and does not need additional differencing.

Now that we know that one difference is necessary, lets move on to a parameter search.

# Define parameter ranges
p_range <- 0:6
d_range <- 1
q_range <- 0:9  

# Calculate total combinations
n_combinations <- length(p_range) * length(d_range) * length(q_range)

# Create an empty matrix to store results
results_matrix <- matrix(NA, nrow = n_combinations, ncol = 6)

# Initialize index for matrix row
i <- 1

# Loop through combinations of ARIMA model parameters
for (d in d_range) {
  for (q in q_range) {
    for (p in p_range) {
  
      # Fit ARIMA model with specified (p,d,q)
      model <- Arima(coneduc_ts, order = c(p, d, q), include.drift = TRUE)
  
      # Store model parameters and AIC/BIC/AICc values in matrix
      results_matrix[i, ] <- c(p, d, q, model$aic, model$bic, model$aicc)
  
      # Increment row index
      i <- i + 1
    }
  }
}

# Convert matrix to data frame
results_df <- as.data.frame(results_matrix)
colnames(results_df) <- c("p", "d", "q", "AIC", "BIC", "AICc")

# Find the row with the lowest AIC
highlight_row <- which.min(results_df$AIC)

# Generate kable table with highlighting for the row with the lowest AIC
knitr::kable(results_df, align = 'c', caption = "Comparison of ARIMA Models") %>%
  kable_styling(full_width = FALSE, position = "center") %>%
  row_spec(highlight_row, bold = TRUE, background = "#FFFF99")  # Highlight row in yellow
Comparison of ARIMA Models
p d q AIC BIC AICc
0 1 0 -146.1814 -142.3177 -145.9314
1 1 0 -154.1756 -148.3801 -153.6650
2 1 0 -160.2939 -152.5666 -159.4243
3 1 0 -161.7735 -152.1144 -160.4402
4 1 0 -160.6396 -149.0486 -158.7305
5 1 0 -159.4146 -145.8918 -156.8099
6 1 0 -161.3818 -145.9272 -157.9532
0 1 1 -157.4889 -151.6935 -156.9783
1 1 1 -155.6315 -147.9042 -154.7619
2 1 1 -164.4317 -154.7726 -163.0984
3 1 1 -162.4820 -150.8911 -160.5729
4 1 1 -160.5157 -146.9929 -157.9110
5 1 1 -159.1041 -143.6495 -155.6755
6 1 1 -159.4157 -142.0293 -155.0255
0 1 2 -155.9474 -148.2201 -155.0778
1 1 2 -158.4396 -148.7804 -157.1062
2 1 2 -162.4920 -150.9010 -160.5829
3 1 2 -162.4265 -148.9037 -159.8218
4 1 2 -160.9624 -145.5078 -157.5338
5 1 2 -159.1660 -141.7795 -154.7757
6 1 2 -159.7447 -140.4264 -154.2447
0 1 3 -158.2586 -148.5995 -156.9253
1 1 3 -159.5933 -148.0023 -157.6842
2 1 3 -162.3082 -148.7854 -159.7035
3 1 3 -161.0628 -145.6082 -157.6342
4 1 3 -158.0646 -140.6782 -153.6744
5 1 3 -156.1666 -136.8484 -150.6666
6 1 3 -159.5780 -138.3280 -152.8088
0 1 4 -162.1637 -150.5727 -160.2546
1 1 4 -162.4235 -148.9007 -159.8188
2 1 4 -160.5475 -145.0929 -157.1189
3 1 4 -159.1997 -141.8133 -154.8095
4 1 4 -158.1940 -138.8757 -152.6940
5 1 4 -156.2203 -134.9702 -149.4510
6 1 4 -160.4167 -137.2348 -152.2062
0 1 5 -162.2605 -148.7377 -159.6559
1 1 5 -161.2509 -145.7963 -157.8224
2 1 5 -159.2674 -141.8809 -154.8771
3 1 5 -159.8130 -140.4947 -154.3130
4 1 5 -158.7977 -137.5476 -152.0284
5 1 5 -157.0933 -133.9114 -148.8828
6 1 5 -155.3629 -130.2491 -145.5250
0 1 6 -160.3046 -144.8500 -156.8760
1 1 6 -159.2570 -141.8705 -154.8667
2 1 6 -157.4700 -138.1517 -151.9700
3 1 6 -156.3756 -135.1255 -149.6063
4 1 6 -157.4528 -134.2709 -149.2422
5 1 6 -156.2534 -131.1397 -146.4156
6 1 6 -157.7621 -130.7165 -146.0954
0 1 7 -160.7412 -143.3547 -156.3509
1 1 7 -158.7809 -139.4626 -153.2809
2 1 7 -160.6976 -139.4475 -153.9284
3 1 7 -158.6980 -135.5161 -150.4875
4 1 7 -156.3615 -131.2478 -146.5237
5 1 7 -154.9472 -127.9016 -143.2805
6 1 7 -155.7651 -126.7877 -142.0508
0 1 8 -158.8138 -139.4956 -153.3138
1 1 8 -159.0816 -137.8316 -152.3124
2 1 8 -158.6981 -135.5162 -150.4875
3 1 8 -156.6980 -131.5842 -146.8601
4 1 8 -157.0020 -129.9564 -145.3353
5 1 8 -152.7827 -123.8053 -139.0684
6 1 8 -154.9958 -124.0866 -138.9958
0 1 9 -161.5029 -140.2528 -154.7337
1 1 9 -159.6613 -136.4794 -151.4508
2 1 9 -158.4304 -133.3166 -148.5925
3 1 9 -152.4592 -125.4136 -140.7925
4 1 9 -155.0022 -126.0249 -141.2880
5 1 9 -153.0854 -122.1762 -137.0854
6 1 9 -153.0750 -120.2339 -134.5295
auto.arima(coneduc_ts)
Series: coneduc_ts 
ARIMA(2,1,1) with drift 

Coefficients:
          ar1      ar2     ma1    drift
      -1.0841  -0.6679  0.6965  -0.0074
s.e.   0.1485   0.1162  0.1436   0.0038

sigma^2 = 0.002027:  log likelihood = 87.22
AIC=-164.43   AICc=-163.1   BIC=-154.77

According to the parameter search and the auto.arima() function, our best model is a simple ARIMA(2,1,1). Lets quickly look at the model diagnostics.

model_output1 <- capture.output(sarima(coneduc_ts,2,1,1))

# Find the line numbers dynamically based on a keyword
start_line <- grep("Coefficients", model_output1)  # Locate where coefficient details start
end_line <- length(model_output1)  # Last line of output

# Print the relevant section automatically
cat(model_output1[start_line:end_line], sep = "\n")
Coefficients: 
         Estimate     SE t.value p.value
ar1       -1.0841 0.1485 -7.3013  0.0000
ar2       -0.6679 0.1162 -5.7458  0.0000
ma1        0.6965 0.1436  4.8502  0.0000
constant  -0.0074 0.0038 -1.9609  0.0558

sigma^2 estimated as 0.001867907 on 47 degrees of freedom 
 
AIC = -3.224151  AICc = -3.207101  BIC = -3.034756 
 

The Residual Plot shows nearly consistent fluctuation around zero, suggesting that the residuals are nearly stationary with a constant mean and finite variance over time.

The Autocorrelation Function (ACF) of the residuals reveals no significant autocorrelations.

The Q-Q Plot indicates that the residuals follow a near-normal distribution, with minor deviations at the tails, which is typical in time series data.

The Ljung-Box Test p-values remain above the 0.05 significance level, implying no significant autocorrelations left in the residuals.

Coefficient Significance: The model coefficients in the ARIMA(2,1,1) model are all significant

Overall, this seems like a strong model. Lets compare it to some benchmark methods.

fit <- Arima(coneduc_ts, order = c(2,1,1), include.drift = FALSE)
summary(fit)
Series: coneduc_ts 
ARIMA(2,1,1) 

Coefficients:
          ar1      ar2     ma1
      -1.0610  -0.6314  0.7172
s.e.   0.1504   0.1203  0.1375

sigma^2 = 0.002134:  log likelihood = 85.43
AIC=-162.85   AICc=-161.98   BIC=-155.13

Training set error measures:
                      ME       RMSE        MAE        MPE     MAPE      MASE
Training set -0.01067867 0.04438487 0.03073701 -0.5315693 1.437563 0.9150038
                    ACF1
Training set -0.01935951
naive <- naive(coneduc_ts,h=5)
mean  <- meanf(coneduc_ts,h=5)
drift <- rwf(coneduc_ts,drift=TRUE,h=5)
accuracy(naive)
                       ME       RMSE        MAE        MPE     MAPE MASE
Training set -0.006753008 0.05591304 0.03359222 -0.3509282 1.552349    1
                   ACF1
Training set -0.4017722
accuracy(mean)
                        ME       RMSE        MAE        MPE     MAPE     MASE
Training set -1.238466e-16 0.09197785 0.06769079 -0.1854871 3.184123 2.015073
                  ACF1
Training set 0.7513554
accuracy(drift)
                       ME       RMSE        MAE         MPE     MAPE      MASE
Training set 3.047601e-17 0.05550374 0.03343843 -0.03171947 1.541556 0.9954219
                   ACF1
Training set -0.4017722

Our ARIMA(2,1,1) beats all the benchmark methods according to the RMSE and the MAE.

Finally, we can plot our forecast.

# Generate the forecast
forecast_result <- forecast(fit, h = 5)

# Plot the forecast
autoplot(forecast_result) +
  labs(title = "ARIMA(2,1,1) Forecast",
       x = "Time",
       y = "Average Confidence in Education") +
  theme_minimal()

The ARIMA predicts that confidence in education will remain low over time.

Conclusion

These analyses imply that that despite the government spending more on education per capita, people believe that education should be prioritized more. One possibility is that government spending on education is being misused and isn’t have as big of an effect on the quality of education as it should. Another possibility is that government spending in general hasn’t been increasing as much as it should so even if the portion of spending that is being used for education increases, the actual amount being spent isn’t meeting the demand.